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Abstract 

This paper is concerned with the Coherent Quantum Linear Quadratic Gaussian (CQLQG) control problem of finding 
a stabilizing measurement-free quantum controller for a quantum plant so as to minimize an infinite-horizon mean square 
performance index for the fully quantum closed-loop system. In comparison with the observation-actuation structure of classical 
controllers, the coherent quantum feedback is less invasive to the quantum dynamics and quantum information. Both the plant 
and the controller are open quantum systems whose dynamic variables satisfy the canonical commutation relations (CCRs) 
of a quantum harmonic oscillator and are governed by linear quantum stochastic differential equations (QSDEs). In order to 
correspond to such oscillators, these QSDEs must satisfy physical realizability (PR) conditions, which are organised as quadratic 
constraints on the controller matrices and reflect the preservation of CCRs in time. The CQLQG problem is a constrained 
optimization problem for the steady-state quantum covariance matrix of the plant-controller system satisfying an algebraic 
Lyapunov equation. We propose a gradient descent algorithm equipped with adaptive stepsize selection for the numerical 
solution of the problem. The algorithm finds a local minimum of the LQG cost over the parameters of the Hamiltonian and 
coupling operators of a stabilizing PR quantum controller, thus taking the PR constraints into account. A convergence analysis 
of the proposed algorithm is presented. A numerical example of a locally optimal CQLQG controller design is provided to 
demonstrate the algorithm performance. 


I. INTRODUCTION 

Coherent quantum feedback control [11], [13] is a relatively novel quantum control paradigm which is aimed at achieving 
given performance specifications for quantum systems, such as internal stability and optimization of a cost functional. Such 
systems arise naturally in quantum physics [7] and its engineering applications (for example, nanotechnology and quantum 
optics [6]). The dynamic variables of quantum systems are (usually noncommuting) operators on an underlying Hilbert 
space which evolve according to the laws of quantum mechanics [15]. The latter make the quantum dynamics particularly 
sensitive to interaction with classical devices in the course of quantum measurement, as reflected in the projection postulate 
of quantum mechanics. In order to overcome this issue, coherent quantum control employs the idea of direct interconnection 
of quantum systems to be controlled (quantum plants) with other quantum systems playing the role of controllers, possibly 
mediated by light fields. Unlike the traditional observation-actuation control loop, this fully quantum measurement-free 
feedback avoids the loss of quantum information as a result of its conversion to classical signals. 

Quantum-optical components, such as optical cavities, beam splitters and phase shifters, make it possible to implement 
coherent quantum feedback governed by linear quantum stochastic differential equations (QSDEs) [17], [18], provided the 
latter are physically realizable (PR) as open quantum harmonic oscillators [5], [6]. The resulting PR conditions [9], [20], 
[21] are organized as quadratic constraints on the coefficients of the QSDEs. The PR constraints for the state-space matrices 
of a coherent quantum controller complicate the solution of quantum counterparts to the classical and Linear Quadratic 
Gaussian (LQG) control problems. 

The Coherent Quantum LQG (CQLQG) control problem [16] seeks for a stabilizing PR quantum controller so as to 
minimize a mean square performance criterion for the fully quantum closed-loop system. A numerical procedure for finding 
suboptimal controllers for this problem was proposed in [16], and algebraic equations for the optimal CQLQG controller 
were obtained in [25]. Despite the previous results, the CQLQG control problem does not lend itself to an “elegant” solution 
(for example, in the form of decoupled Riccati equations as in the classical case [10]) and remains a subject of research. 
Since the main difficulties are caused by the coupling of the equations due to the PR constraints, a conversion of the CQLQG 
control problem to an unconstrained problem by using Lagrange multipliers was considered in [26] for a related coherent 
quantum filtering problem which is a simplified feedback-free version of the CQLQG control problem. 

In the present paper, we develop an algorithm for the numerical solution of the CQLQG control problem by using the 
gradient descent method and the Hamiltonian parameterization of PR quantum controllers [25]. The latter is a different 
technique to handle the PR constraints by reformulating the CQLQG control problem in an unconstrained fashion. More 
precisely, the optimal solution is sought in the class of stabilizing PR controllers whose state-space matrices are parameterized 
in terms of the free Hamiltonian and coupling operators of an open quantum harmonic oscillator [5]. We obtain ordinary 
differential equations (ODEs) for the gradient descent in the Hilbert space of matrix-valued parameters. Lor this purpose, 
Lrechet differentiation is used together with related algebraic techniques [3], [14], [22], [23], [25] to employ the analytic 
structure of the LQG cost as a composite function of the matrix-valued variables, involving Lyapunov equations. The 
advantage of the proposed approach is that, at intermediate steps, it produces stabilizing PR quantum controllers which can 
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be used as gradually improving suboptimal solutions of the CQLQG control problem, and a locally optimal solution (if it 
exists) is achieved asymptotically by moving along negative gradient directions with a suitable choice of stepsizes. To this 
end, we provide an algorithm for adaptive selection of the stepsize for each iteration based on the second-order Gateaux 
derivative of the LQG cost along the gradient. However, the proposed gradient descent algorithm for the CQLQG control 
problem requires for its initialization a stabilizing PR quantum controller. Finding such a controller for an arbitrary given 
quantum plant is a nontrivial open problem. Because of the lack of a systematic solution for this quantum stabilization 
problem at the moment, the current version of the algorithm is initialized at a stabilizing PR quantum controller obtained 
by random search in the space defined by the Hamiltonian parameterization of PR controllers. Although a random search 
for an admissible starting point is acceptable for low-dimensional problems, the development of a more systematic solution 
for this issue is a subject of future research. 

The paper is organised as follows. Sectionj^outlines the notation used in the paper. Sections [^andj^specify the quantum 
plants and coherent quantum controllers being considered. Section |V] revisits PR conditions for linear quantum systems. 
Section [Vl|formulates the CQLQG control problem. Section [Vn] describes a gradient descent system for finding local minima 


in the control problem. Section VIII describes an algorithmic implementation of the gradient descent method. Section IX 


discusses convergence of the algorithm. Section provides a numerical example of designing a locally optimal CQLQG 
controller. Section XI gives concluding remarks. Appendices [A| and [B] provide a subsidiary material on the differentiation 
of the LQG cost. 


11. NOTATION 

Vectors are assumed to be organized as columns unless specified otherwise, and the transpose acts on matrices with 
operator-valued entries as if the latter were scalars. For a vector X of operators Xi,...,Xr and a vector Y of operators 
Fi,...,ys, the commutator matrix [X,y^] :=XY^ —{YX^Y (r x s)-matrix whose {j,k)th entry is the commutator 

[Xy,y)-] '■= ^jYk — YkXj of the operators Xj and Yk- Furthermore, := ((•)*)^ denotes the transpose of the entry-wise 
operator adjoint {■)*■ When it is applied to complex matrices, reduces to the complex conjugate transpose (•)* := ((O)^- 
Denoted by sym(-) := and asym(-) := are the symmetrizer and antisymmetrizer of matrices. Also, we denote 

by Sr, Ar and Hr := Sr-f/Ar the subspaces of real symmetric, real antisymmetric and complex Hermitian matrices of 
order r, respectively, with i := v^—T the imaginary unit. Furthermore, Ir denotes the identity matrix of order r, positive 
(semi-) definiteness of matrices is denoted by (^) and 0 is the tensor product of spaces or operators (in particular, the 
Kronecker product of matrices). The adjoints and self-adjointness of linear operators acting on matrices is understood in 
the sense of the Frobenius inner product {M,N) : = Tr{M*N) of real or complex matrices, with the corresponding Frobenius 
norm ||M|| := Y{M,M) which is the standard Euclidean norm | • | for vectors. Also, EX := Tr(pX) denotes the quantum 
expectation of a quantum variable X (or a vector of such variables) over a density operator p which specifies the underlying 
quantum state. 

III. QUANTUM PLANT 

The quantum plant under consideration is an open quantum stochastic system which is coupled to another such system 
(playing the role of a controller), with the dynamics of both systems being affected by the environment. Both the plant and 
the controller are assumed to satisfy the physical realizability (PR) conditions [9], [16], [20] which will be described in 
Section|v] The plant has n dynamic variables xi(f),... (with n even) which are self-adjoint operators on a Hilbert space 
specified below. With the time arguments being omitted for brevity, the evolution of the plant state vector x := {xk)i^k^n 
and its contribution to a pi-dimensional output of the plant y := {yk)i^k^pi (also with self-adjoint operator-valued entries) 
are governed by QSDEs 

dx = Axdt + Bdw + Edi], dy = Cxdt-b Ddw. (1) 

Here, A G B G C G D G E G are given constant matrices. Also, 

z := Cx (2) 

is a “signal part” of the plant output y, and Tj is a p 2 -dimensional output of the controller to be described in Section 
The external noise acting on the plant is represented by a quantum Wiener process w := (wk)i^k^mi whose entries are 
self-adjoint operators on a boson Eock space [17] with the quantum Ito table dwdw^ = fiidf, where the matrix fii G H„,| 
is given by fii 7^, +/7i ^ 0. Here, the matrix 7i G A,„^ specifies the CCRs between the entries of the plant noise w as 

[dw,dw^] = 2iJ\dt and (assuming that the dimension m\ is even) is given by J\ := /m ,/2 ® X with J ;= 


0 1 

-1 0 
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IV. QUANTUM CONTROLLER 

Consider an interconnection of the plant ([T]) with a coherent (that is, measurement-free) quantum controller. The latter is 
also an open quantum system with an n-dimensional state vector ^ := {^k)\^k^n of self-adjoint operators on another Hilbert 
space, which also evolve in time. The assumption that the controller has the same number of dynamic variables as the plant 
is adopted from the classical LQG control theory. The controller interacts with the plant and the environment according to 
the QSDEs 

di* = a^df-f hdo)-fedy, drj = c^dr-l-t/da). (3) 

Here, a € b € ^ g j g g g K«xpi constant matrices. Similarly to (|^, the p 2 -dimensional 

process 

C-=c^ (4) 


is the signal part of the controller output rj. The process O) in Q is a quantum noise which effects the controller and is 
an m 2 -dimensional quantum Wiener process (with m 2 even) on another boson Eock space ^2 with the quantum Ito table 
dodo)^ = f22dt, where the matrix Q .2 G is given by Q .2 := /mj +iJ 2 ^ 0. Here, the matrix J 2 G specifies the CCRs 
between the entries of the controller noise (O as [da),da)^] = 2 iJ 2 dt and is given by J 2 := plant and controller 

noises w and a> act on different boson Eock spaces and ^27 respectively, and hence, commute with each other. Therefore, 
the combined quantum Wiener process 

w 


W:= 


0 } 


of dimension m\=m\+m 2 acts on the tensor product space ®,^2 and has a block diagonal CCR matrix J: 


[dW,AW'^]=2iJdt, 


J- 


Ji O' 
0 J 2 


(5) 


(6) 


Eurthermore, the external boson fields are assumed to be in the product vacuum state v := Vi0V27 and hence, are uncorrelated. 
The resulting quantum Ito table of the combined Wiener process W'm ^ '\s 

= Q.dt, Q.:=I,r, + iJ = n* > 0. (7) 


In the controller dynamics ([^, the matrix b is the noise gain matrix, while e plays the role of the observation gain matrix, 
although y is not an observation signal in the classical control theoretic sense. Accordingly, the process 1 ^ in (j^ corresponds 
to the classical actuator signal. Now, the combined set of QSDEs (jT]), Q describes the fully quantum closed-loop system 
shown in Eig. [T| By using a quadratic cost adopted in quantum control settings [16], [25] from classical LQG control [10], 
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Fig. 1. The interconnected quantum plant and quantum controller form a fully quantum closed-loop system which is governed by (D'S and is influenced 
by the environment through the quantum Wiener processes w, ft). 


the performance of the coherent quantum controller will be described in Section VI in terms of an r-dimensional quantum 
process 

3r:=Fx + GC. (8) 


Here, F G G G are given weighting matrices whose entries quantify the relative importance of the state variables 

xi,...,Xn of the plant and the “actuator output” variables ^ 1 ,..., of the controller. The choice of F, G is specified only 
by control design preferences and is not subjected to physical constraints. The process .S° in ([^ is linearly related to the 
2n-dimensional state vector ^ , 


jr:= 


X 


(9) 


of the closed-loop system whose dynamics are governed by the QSDE 


d^ = £/^dt + mW, = 


( 10 ) 
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which is driven by the combined quantum Wiener process W from (jsjl. The state-space matrices 

g Qf (}jg closed-loop system in (10 1 are obtained by combining with (j^, (j^ and depend on the controller 

matrices a, b, c, e in an affine fashion: 


: = 


A 

eC 


Ec 

a 


B 

eD 


Ed 

b 


'■.= [F Gc\. 


( 11 ) 


V. CONDITIONS FOR PHYSICAL REALIZABILITY 

Both the quantum plant ([T]| and the coherent quantum controller Q are assumed to be physically realizable as open 
quantum harmonic oscillators, with initial complex separable Hilbert spaces In particular, their dynamic variables 

(which are self-adjoint operators on the product space at any subsequent moment of time f > 0) satisfy 

CCRs 

[x,x^]=2i®u K,^^]=2i02, (12) 

where ©i,02 € A„ are constant nonsingular matrices. An equivalent form of the CCRs for the combined vector ^ from 
(|^ is 


[S', S''^] = 2/0, 


0 := 


01 

0 


0 

02 


(13) 


The preservation of the CCRs (12 1 (including the commutativity between x and ^) is a consequence of the unitary evolution 
of the isolated system formed from the plant, controller and their environment. The QSDE in ([T0|) preserves the CCR matrix 
0 in (13 I in time if and only if the matrices S in (111 satisfy 




(14) 


where J is the CCR matrix of the combined quantum Wiener process W in 0 given by (j^. The relation ( [l4l l is obtained 
by taking the imaginary part of the algebraic Lyapunov equation (ALE) 

S^S + S£/'^ + = 0 (15) 

(provided £/ is Hurwitz) for the steady-state quantum covariance matrix 

5:=^lim E(S'(f)S'(f)’^) = P-f/0 = 5* ^ 0, (16) 

with Q. the quantum Ito matrix from 0. Here, the quantum expectation E( ) is taken over the product state (U0V, where 
Iff is the initial quantum state of the plant and controller on 0 ,Jf 2 , and v is the vacuum state of the external fields on 
,^10^2- We have also used the convergence lim(^+„oE,^(f) =0 which is ensured by being Hurwitz. The real part 

F:=ReS (17) 

of the quantum covariance matrix S from is the unique solution to the ALE 

+ + = 0, (18) 

obtained by taking the real part of ([Ts]), and coincides with the controllability Gramian [10] of the pair (.c/,^). Since 


(19) 

( 20 ) 
( 21 ) 


the left-hand side of (14 1 is an antisymmetric matrix of order 2 n, then, by computing the diagonal (n x n)-blocks and the 
upper off-diagonal block of this matrix with the aid of ( |TT] i, it follows that the preservation of the CCR matrix 0 in ( |T3] l is 
equivalent to 

A&i+&iA'^ + BJiB'^+EdJ2d'^E'^ = 0, 

a&2 + &2a^ + eDJiD^e^ + bJ2b^ = 0 , 

{@iC^ +BJiD^)e^ + E{c@2 + dJ2b'^) = 0\ 


cf. [24, Eqs. (18)-(20)]. Therefore, the fulfillment of the equalities 

0iC'^-t-B7iD’^ = O, 
C02 -f dJ2i^ = 0 


( 22 ) 

(23) 


is sufficient for 0. Note that ( [T9] l, are the conditions for physical realizability (PR) [9], [16], [20] of t he q uantum 
plant which describe the preservation of the CCR matrix 0i in (12 1 and [x,y^] = 0. Similarly, the relations (20 1 , (23 i, which 
describe the preservation of the CCR matrix 02 in (12 1 and = 0, are the PR conditions for the coherent quantum 

controller. The PR condition (20 1 can be regarded as a linear equation with respect to the matrix a, and its general solution 
is representable as 

^ -iT T I 7,1 LT'ir',-! 


1 


a = 2@2R--{eDJiD^e' +bJ2b')&2 


(24) 
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Here, the matrix /? € S„ specifies the free Hamiltonian which the PR controller would have in the absence of 

interaction with its surroundings; cf. [5, Eqs. (20)-(22) on pp. 8-9]. The other PR condition ( [2^ allows the matrix c to be 
expressed in terms of b as 

c = -dJ2b'^®2''- ( 25 ) 

The coupling between the output matrix c and the noise gain matrix b makes the design of a coherent quantum controller ([^ 
substantially different from that of the classical controllers even at the level of achieving internal stability for the closed-loop 
system. Indeed, if the additional quantum noise (O is effectively eliminated from the state dynamics of the quantum controller 
by letting b — Q, then 


implies that c = 0 , and hence, the matrix £/ in ( 111 becomes block lower triangular. In this case, 
the closed-loop system in ( [T0| ) cannot be internally stable if A is not Hurwitz. Also note that, in the formulations of the PR 
conditions [9], [16], [21] for the plant and controller QSDEs Q, Q, the noise feedthrough matrices are usually given by 
D = [/pj O] and d = [ip^ O], with pi ^ nii and p 2 ^ m 2 . Such matrices D and d have full row rank and satisfy 


DD^ = In 


dd^ = In 


(26) 


Pj, uu — ip2- 

The full row rank property of D corresponds to nondegeneracy of measurements in the classical setting, where y in ([T]) is 
an observation process. Eurthermore, since det 72 7 ^ 0 and det 02 7 ^ 0, the full row rank property of d implies that the map 
9 I— c € given by (25 1 , is onto. This allows the matrix c to be assigned any value by an appropriate choice 

of b, which plays a part in the stabilization issue mentioned above. Although (j2^ simplifies the algebraic manipulations, it 
is the rank properties of the matrices D, d that are most important. 

VI. COHERENT QUANTUM LQG CONTROL PROBLEM 
Eollowing [16], [25], we formulate the CQLQG control problem as that of minimizing the steady-state mean square value 


:= i^lim E(ir(f)Tir(f)) = 


(27) 


for the criterion process ^ of the closed-loop system ( [TOl i over internally stabilizing (that is, making the matrix Hurwitz) 
PR quantum controllers ([^ of fixed dimensions described in Sections IV 10 Here, = YJk=\ squared 

entries of iF (and hence, is a positive semi-definite self-adjoint operator) and P is the controllability Gramian of the 

closed-loop system given by ( [T8] i. The LQG cost S in ( [27] i is a function of the triple 


u := {R,b,e) G S„ x x =: U (28) 

which parameterizes PR quantum controllers (j^ through (24i, (25 1 , with the controller noise feedthrough matrix d G IRP 2^'"2 
being fixed and satisfying (p^. Accordingly, the minimization in (p7jl is carried out over the set 


Uo := {m G U : . 0 ?' in (11) is Hurwitz} 


(29) 


of those u which specify internally stabilizing PR quantum controllers for the quantum plant ([T]). Eor what follows, the 
set U on the right-hand side of ( [28] l is endowed with the structure of a Hilbert space with the direct sum inner product 
{{R,b,e),{R',b',e')) := {R,R') + {b,b') + {e,e’). Note that Uq in (29 1 is an open subset of U. 


VII. GRADIENT ELOW EOR THE LQG COST 

The gradient descent approach to the solution of the CQLQG control problem is to move with the negative gradient flow 
for the LQG cost function S’ in towards a local minimum. The gradient descent can be regarded as a dynamical system 
governed by the ODE 

u{T:) = -8iuiT)), m(0) = mo. (30) 

Here, () := is the derivative with respect to fictitious time T ^ 0, and the gradient 

g{u) := duS{u) = {dRS,dbS,deS) (31) 

is the Erechet derivative of the LQG cost with respect to u in the Hilbert space U associated with the Hamiltonian 


parameterization of PR quantum controllers in (28 1 . More precisely, the map g : Uo —7 U is well-defined on the open 


set Uo in (29 1 . The starting point in (30i is assumed to satisfy 

Mo := (Ro,^o,eo) G Uo, 


(32) 


so that the corresponding PR controller is internally stabilizing. Unless uq is a stationary point of S, the LQG cost is strictly 
decreasing along the trajectory of the ODE (30 1 in view of S{u{z)) = —|jg(M(T))|p. The first-order Erechet derivative in 


(31 1 is computed in the following lemma. Eor its formulation, we denote by Q the observability Gramian of the pair (x/, 'S) 

+ + = (33) 


which is a unique solution of the ALE 
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provided the matrix in (11 1 is Hurwitz. Furthermore, we will use the Hankelian of the closed-loop system defined by 


H-.= QP. 


(34) 


Also, we partition (2n x 2n)-matrices X (such as P, Q, H) into (n x n)-blocks Xjj^ as X := 
Lemma 1: For any m G Uq from 


Xl 

^ 12 ' 

t 

A 21 

^22. 



, the Frechet derivative (j3T]l of the LQG cost S' m {21 \ can be computed as 

SrS = -2sym(02//22), 

= QixEd + Q22b - ij/bJi- XdJi, 
dgS = H2iC^ -f Q2iBD^ + 022^ ~ \j/eDJiD^. 

Here, \j/ and x auxiliary (n x n)-matrices defined by 

V/ ;= asym(//22©2 '), I := ©2 ’ + ^’ 21 ^’'^^ + P22c'^G^G), 


(35) 

(36) 

(37) 

(38) 

with P, Q, H the Gramians and Hankelian from ( [T8] l, ( |3^ , ( |5^ . □ 

The proof of Lemma is similar to that of [25, Theorem 1] and is given in Appendix [A| for completeness. That the 
trajectories of the gradient descent system in ( [30l l will not “miss” local minima of the LQG cost is justified by the following 
lemma. 


Lemma 2: A point m* G Uq in (29 i is a stable equilibrium of the ODE (30 1 if and only if it is a local minimum of the 
LQG cost S in □ 

Proof: The assertion of the lemma can be established by using [1, Theorem 3] and the analyticity [ 8 ] (rather than 
inhnite differentiability) of the LQG cost S in a neighbourhood of any point u G Uq. The analyticity follows from the 
representation 


S = — (£/ Q£/) 


vec 




(39) 


which is obtained from (I 81 , (271 by using the column-wise vectorization vec(-) of matrices [14], [22] and the Kronecker 
sum sS := l 2 n ^ s/ + ^ l 2 n of the matrix j:/ with itself. Indeed, the representation (39i implies that ^ is a rational 


function of the entries of £/, ^ in (11 1 which, in turn, are polynomial functions of the entries of R, b, e in view of 


(24 1 , (251, and hence, S{u) is a rational function of u. Therefore, the function S{u) is analytic on the open set Uo since the 


matrix 22 / 0 . 2 / is also Hurwitz (and hence, nonsingular) for any u G Uq. ■ 

In practice, the gradient descent ODE ( [30| is solved by using a numerical algorithm, which is the subject of the next 
section. 


VIII. GRADIENT DESCENT ALGORITHM 

We will now consider a numerical algorithm which implements the gradient descent method ( [30| ) for the CQLQG control 
problem in the form 

Uk+i-=Uk-Sk8iuk), k = Q,l,2,.... (40) 


This recurrence equation is initialized with matrices Rq, bo, eo of an internally stabilizing PR controller in (32i (see 
VIII-AI. The gradient g{uk) is computed by using Lemma [T| and the stepsize > 0 is chosen as described in 
VIII-BI The iterations in (|40|) are stopped when a termination condition is satisfied (see Section |VIII-C|l. The 


Section 
Section 

ingredients of the algorithm are discussed in the subsequent sections. 
A. Initialization 


The gradient descent algorithm (40i relies on existence of an internally stabilizing PR quantum controller which can be 
used as an initial point. As mentioned in Introduction, the existence of such controllers (that is, nonemptiness of the set 


Uq in (29 1 ) for a given quantum plant (and a systematic method of finding them) remains an open problem. In the present 
version of the algorithm, this quantum stabilization problem is solved by using a random search in the finite-dimensional 
Hilbert space U in 
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B. Stepsize selection 

According to the conventional limited minimization rule (see, for example, [4]), the stepsize Sk is chosen for each iteration 
of the gradient descent by solving the minimization problem 


Sk G AYgvamS{uk-sg{uk)) 

O^s^hk 


(41) 


with a constant search horizon hk := h> 0. Here, we use the convention that ^(m) := +oo if m 7 ^ Uq (thus discarding those 
controllers which are not internally stabilizing). A restricted version of the line search with a constant horizon h may suffer 
from the inability to adapt properly to the behaviour of the function S’ in its minimization over the ray {uk — sg{uk) : s > 
0} C U. In order to overcome this issue, for the stepsize selection in the gradient descent algorithm ( |4^ , we will use a 
modihed version of the limited minimization rule with an adaptive choice of the search horizon hk in each iteration. More 
precisely, hk can be chosen so as to enable (411 to “capture” the minimum of S over the whole ray if tl’ is a strictly convex 


quadratic function. To this end, consider the quadratically truncated Taylor series 


where 


S{u — sg) = S[u) — s^gS + —^gS + o(s^), 


%S{u) := + = {g{u),v), 

^^S{u) := d}S{u + sv)\ = {dlS{u){v),v) 


(42) 

(43) 

(44) 


are the hrst and second-order Gateaux (or directional) derivatives [12] of the LQG cost at a point m G Uo (specifying an 
internally stabilizing controller) along v G U. Here, in view of ( |43| ), 

%S = \\g\\^ = \\dnS\\^ + \\dkSf + \\d,S\\'^- 


^ 0 . 


(45) 


Also, d^Sju) := dug(u) in ( |44| is the second-order Frechet derivative of S which is a self-adjoint operator on the Hilbert 
space U in whose computation is outlined in Appendix [ b| Now, if StgS{u) > 0, then the quadratic polynomial of s 
on the right-hand side of (42i (with the higher-order terms being neglected) achieves its unique minimum at a nonnegative 
value of s: 


argmin 




Ikll^ 


SijS' 


(46) 


This suggests using the right-hand side of (46 1 as a search horizon hk in (411, provided S!gS(u) > 0. However, if the latter 


inequality does not hold, the argument, basea on a quadratic approximation of the minimization problem (41 1 , is no longer 
valid and needs to be amended. In this case (when @'^S{u) ^ 0), the search horizon can be chosen so as to avoid the 
domination of nonlinear terms over the linear term in the quadratically truncated Taylor series for the LQG cost along the 


ideal gradient descent trajectory in (30i: 


S{u{t + s)) = S(u{t)) + {S{u)) s + {S{u)) 


= S{u{T))-\\g\\h- 


-^ISs^- 




(47) 


where (44i is used. For i ^ 0, the comparison of the absolute values ||g|pi and \Sl'^S\s'^ of the linear and quadratic terms 
in ( |47| ) shows that the latter does not dominate the former if 

II 

5^741^4:7. (48) 


\9lS\ ■ 

This inequality is closely related to the accuracy of ((40|l as Euler scheme for numerical integration of the ODE ([30|. More 
precisely, if the stepsizes > 0 in (401 are signihcantly smaller than the respective values of the right-hand side of (48 1 , 
then Uk becomes an accurate approximation of the ideal gradient descent trajectory m(t) at hctitious time T ; = io -f.. . + Sk-i. 
A combination of (46i and ( [48| justihes the following heuristic rule for choosing the search horizon at the current point 
Uk G Uq: 


hk := min h^ 


kl4 


I ^2^1 


(49) 


U=Ui( 


Here, h„ 


1 a given positive threshold which becomes active, for example, if ^^S vanishes. The stepsize selection algorithm 
considered below, replaces the minimization problem in ( [4T] i with a different procedure which involves a hnite subset of 
values of s from a geometric progression 

Sk/--=hkf, £ = 0,1,2,... (50) 
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whose initial value is given by (49 1 . The common ratio / € (0,1) is a parameter of the algorithm which affects how 
“densely” the progression fills the interval [0,/!^]. Now, the adaptive stepsize selection algorithm proceeds as follows: 

•— ^k. ji ( 51 ) 


where the y'th element of the geometric progression (50i is chosen according to the Armijo rule [4] with a parameter 
(Te(0,l): 


;■ :=min{^^0: S{uk) - S{uk-Sk/g{uk)) ^ ask/\\g{uk)f}. 


(52) 


Here, the subset of indices i is nonempty since ff < 1 and lim^_^+„ = 0. The inequality in (521 is important in the 

convergence analysis of the gradient descent algorithm. In particular, the condition a > 0 ensures that S‘{uk) is strictly 
decreasing until achieves a stationary point of the LQG cost. Such a point is a stable equilibrium of the gradient descent 
only if it delivers a local minimum to the LQG cost. 

C. Termination condition 


Since the gradient descent sequence in (401 can converge to a local minimum of the LQG cost only asymptotically, as 
k —> + 00 , the algorithm is equipped with a termination condition (for stopping the iterations) which reflects the proximity 
to the limit point. More precisely, we use the following termination condition which employs the relative smallness of the 
gradient as specified by a dimensionless parameter £ > 0: 

■S<:|k(Mi:)IK e||«i:||- (53) 

IX. CONVERGENCE OE THE GRADIENT DESCENT ALGORITHM 
As the proposed algorithm is based on the classical gradient descent approach, its convergence analysis follows a similar 
reasoning, which we provide below for completeness. 


Theorem 1: Suppose {uk)k^o is the gradient descent sequence in (40l with the stepsize selection described by (49i-(52i. 
Then every limit point m* G Uq of this sequence is a stationary point of the LQG cost that is, gjw*) = 0. 

Proof: Since the sequence S'{uk) 0 is nonincreasing, it has a finite limit. Therefore, (40i, and the Armijo rule 
in (52i imply that ask\\g{uk)\\'^ ^ ^{uk) — ^{uk+i) —>■ 0 as A: —)■ +oo. Hence, in view of O’ > 0, it follows that 


= 0 . 


(54) 


Now, assume that the gradient descent sequence has a limit point m* := limj(r^k^+ooUk G Uq such that g(ut) f 0, where 
JP := {Q ^ ki < k 2 < ...} is an infinite subset of nonnegative integers which specifies the respective subsequence of 
Then the analyticity of the LQG cost on the open set Uo implies that 


lim g{uk)=g{u^)fQ, 
lim hk = min /zmax 

JG3k^+^ \ 


l|g(«*)lP 


> 0 , 


(55) 

(56) 


where use is made of (49 1 . Note that, if = 0, the limit in (561 is equal to /zmax > 0. A combination of (55i with 

([54|l implies that 


lim Sk = 0. 

JG3k^+oc 


(57) 


In turn, by combining ( [57] i with ( [56] l and recalling 0 and the condition 0 < / < 1, it follows that the indices jp := 
log^ = — '‘linf of elements of the geometric progression in (501, which correspond to kp G diverge to infinity 




2 lS p ^ + 00 , and hence, jp ^ 1 for all sufficiently large p. For all such p, the stepsize candidates Sk^j^-i = -^ do not pass 
the Armijo selection rule in (521, that is. 


S’{uk)-S'{uk-jg{uk)) < aj\\g{uk)f 


Jk, 


(58) 


for all sufficiently large kG JP. Upon multiplying both parts of (58 1 by ^ and taking the limit, this inequality leads to 


llg(“*)ll^ = 

^ (J lim 


(^y^(p'iuk)-p'{uk-jg{uk)) 


\\giuk)\\ = (y\\g{u* 


(59) 


where use is also made of (55i and (57i. However, since O’ < 1, the inequality in (59 1 contradicts the assumption that 
g{ut) 7 ^ 0. This contradiction shows that any limit point m* G Uq of the gradient descent sequence satisfies g(ut) =0. ■ 
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Note that the CQLQG control problem inherits a special type of symmetry from the LQG cost S’ which is invariant under 
symplectic similarity transformations of the controller variables ^ i—> with E G satisfying E02E^ = 02 and thus 

preserving the CCR matrix 02 (see, for example, [24], [25]). Hence, the stationary points of the LQG cost are not isolated, 
which complicates the convergence rate analysis for the proposed gradient descent algorithm. This issue is beyond the scope 
of the present study and will be addressed elsewhere by using more advanced analytic tools (such as in [2] and related 
references). 


X. A NUMERICAL EXAMPLE OE OPTIMAL CQLQG CONTROLLER DESIGN 


The gradient descent algorithm of Section VIII was tested to find a locally optimal solution of the CQLQG control problem 
for a PR quantum plant in (jT]) with dimensions n = m 2 = p\ =p 2 = r = 2,m\ =4 and randomly generated state-space matrices 
A, B, C, E, satisfying the PR conditions ([T9]l, (|22)i, and the weighting matrices F, G in 


0.9534 

-1.1165 

, B = 

■-1.7174 

-0.2189 

1.9180 

0.5636 ■ 

0.4193 

1.8821 

-0.6815 

1.3570 

0.2985 

-0.3679 


C = 

■-1.3570 
-0.6815 

-0.2189 

1.7174 


1 0 0 0' 
0 10 0 

5 

E = 

'-0.3238 0.2779 ' 

-1.1693 -0.5966 

F = 

-0.8290 

-1.8655 

-0.9665' 

-0.0357 

, G = 

■-0.2324 - 
-0.5822 - 

0.1608' 

1.0961 

, d = 

1 0' 
0 1 



This plant is unstable (the eigenvalues of the matrix A are 1.4177±0.5025/). The algorithm was run with parameters /imax = 1, 
/ = 0.5, (7 = 0.9, e = 10^® in (491-(53 i for 10 randomly generated stabilizing PR controllers as initial points. Starting from 
these points, it has taken 307 to 2318 steps for the algorithm to reach the fulfillment of the termination condition, with the 
average number of iterations being 1075. The local minimum value of the LQG cost is #niin = 12.1026 and is achieved at 
the following controller parameters: 


R = 

'-0.5611 

-1.5567' 


' 1.8111 

0.7201 ' 


'-0.1250 

4.9673 ■ 

-1.5567 

1.8283 

, b = 

-1.4979 

-3.9696 


-4.4929 

-1.3387 


The values of the LQG cost S{uk) for the gradient descent sequences are presented in Pig. in the form of semi- 
logarithmic graphs of — 1. These graphs are in qualitative agreement with the relatively slow linear convergence rate, 

®niin 



Fig. 2. The relative deviations — 1 from the minimum value of the LQG cost on a logarithmic scale versus the number of steps k. 

typical for gradient descent methods. However, they also show that the proposed algorithm is fairly reliable, being able to 
cope with poor initial approximations where the LQG cost exceeds the minimum value by an order of magnitude. 

XL CONCLUSION 

A gradient descent algorithm has been developed for the numerical solution of the optimal CQLQG controller design 
problem, and its convergence has been investigated. The algorithm has been tested and it appears to be fairly reliable in 
a numerical example with randomly generated stabilizing PR controllers as initial points. The lack of a more systematic 
method for initialization and a relatively slow convergence rate are the main shortcomings of the algorithm. These issues 
are a subject for future research and will be tackled in subsequent publications. 
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Appendix 

A. Computing the Frechet derivative of the LQG cost 

From the ALEs ( [3^ and the properties of the Frobenius inner product of matrices, it follows that the LQG cost ^ 
in (27 1 is representable as 






(Al) 


where H is the Hankelian given by ( [3^ . In what follows, 5( ) denotes the first variation, and 5x( ) is the first variation with 
respect to an independent matrix-valued variable X. Since the matrix R influences the LQG cost S' in ( |Al| l only through the 
controller matrix a in ( |24l i which enters the matrix sS in ( [TT] i, then 

5rsJ- 5rP) = -^-{sS'^Q + QsS, 5rP) 


= - {Q,sS5RP+i5RP)sSy = -{Q,{5 r£/)P + P5rsS 


2 

={H,Sr 


= (7722,5/ja) = (7722,20257?). 


(A2) 


Now, since /? S S„, and the subspaces S„ and A„ are orthogonal, then the first variation of the LQG cost in (A2i takes the 
form 5rS = 2(027/22,57?) = —2(sym(027/22),57?), which, by the definition of the Frechet derivative, establishes (351. By 
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a similar reasoning, the first variations of the LQG cost with respect to b and e are as follows: 


5b<g ={H, 5b£^) + {Q.SS, dS§) + 5b^) 

= {H22,Sba) + {E^Hi 2, dbc) + {{Q^)22, Sb) 

+ {G^FPn + G^GcP22,5bc) 

= ( - asym{H22&2^)bJ2 + Q2iEd + 022^ 

- Q2\Hj2E+P2iF'^G + P22c'^G^G)dJ2,5b), 

5e(o ={H2\,5e C) + {H22,5ea) + -{Q,5e{^^^^)) 

= {H 2 iC^, 5e) + (// 22 © ',asym(5eD7iD^e^)) 

+ {Ql^,5e£§) 

= (//2iC’^ - asym(//22©^')e£)7iZ)T 

+ {Q2iB + Q22eD)D^,de), 

which leads to ( [36l l and ( |J7| ) in view of ( |26l l and ( |38| ). 


B. Computing the second-order Gateaux derivative of the LQG cost 

The first-order Gateaux derivative of the LQG cost S' along the gradient g in (31 1 is expressed in terms of the first-order 
Frechet derivatives from (35 l-([J7| by using ( [45] l, provided u G Uq. Hence, the second-order Gateaux derivative along the 
gradient can be computed as 


^lS = %{\\g\?)=^{%8,g) 

=2{{^gdRS,dRS) + {^gdbS,dbS) -P {^gdeS,deS)), (Bl) 


where 


9gdRS =-2sym(02^g//22), (B2) 

9gdbS =9gQ2xEd + 9g{Q22b) 

-&g{\l/b)J2-9gXdJ2, (B3) 

^gdeS =9gH2xC'^ + ^gQ2\BD^ 

+ ^g{Q22e)-%{we)DhD^. (B4) 


The second-order Gateaux derivative in •ET} can now be computed by using ( |B2|l-(|B4| l, the Leibniz product rule and the 
first-order Gateaux derivatives of P, Q. By differentiating both sides of the ALE (18i and its dual ( |3^ , it follows that the 
'’gQ are unique solutions of the ALEs 


matrices ^gP and 


sSS>gP -\- ^gPsi'^ -I- l'sym{SigsSP S = 0, 

sS'^^gQ -\- !^gQsS -\- 2^ym{sS^^gQ -\- = 0. 


Here, in view of (111, (24 1 , (251 and the relation ^gU = g. 


^gjS = 




0 -EdJ2dbST'&2^ 
dgSC Siga 

■ 0 o' 

dgSD dbS ’ 


^g'S=[0 -GdJ2dbS'^®f'^], 


with 

3)ga = 202(9«(f — asym((9(,(^’L>/iL>^e^ -(-<3/,<^/2^^)©2'• 
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